#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jul  6 12:11:08 2023

@author: jcfq2
"""

"""Now let's flip to System III W coordinates"""
import JupiterMag as jm 
import matplotlib.colors as colors
import matplotlib.pyplot as plt

from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np










def Bmag(longc,latgc,model='jrm33',radial=True,R=1.00,MaxDeg=30):

    
    longc2 = 360-longc
    latgc2 = latgc+90
    
    
    longcr = longc2*np.pi/180.0
    latgcr = (180-latgc2)*np.pi/180.0
    r = np.zeros(longcr.shape) + R
    
    jm.Con2020.Config(equation_type='analytic',CartesianIn=False,CartesianOut=False)
    jm.Internal.Config(Model=model,CartesianIn=False,CartesianOut=False,Degree=MaxDeg)
    BrI,BtI,BpI = jm.Internal.Field(r,latgcr,longcr,MaxDeg=MaxDeg)
    BrE,BtE,BpE = jm.Con2020.Field(r,latgcr,longcr)
        
    B = np.sqrt(BrI**2 + BtI**2 + BpI**2) + np.sqrt(BrE**2 + BtE**2 + BpE**2)
    
    #convert to Gauss
    Br = BrI + BrE
    if radial:
        Bg = Br.reshape(longcr.shape)*1e-5
    else:
        Bg = B.reshape(longcr.shape)*1e-5

    return Bg



latt = np.linspace(180,0,181)
long = np.linspace(0,360,361)
latc = 0.5*(latt[1:] + latt[:-1])
lonc = 0.5*(long[1:] + long[:-1])
longc,latgc = np.meshgrid(lonc,latc)

latitude = 180-(latgc)
longitude = longc


Bg = Bmag(longc,latgc,model='jrm33',radial=True)


#get the scale
scale = [-60.0,60.0]
    
#set norm
norm = colors.Normalize(vmin=scale[0],vmax=scale[1])
        
maps = [1,1,0,0]

fig = plt.figure(figsize=(15,              12))
ax = plt.subplot(111)
ax.set_aspect(1.0)
ax.set_xlabel('SIII West Longitude ($^\circ$)',              fontsize=14)
ax.set_ylabel('SIII Latitude ($^\circ$)',              fontsize=14)
        
plt.title('Iso-contours of the strength of the radial component of the field (JRM33 + Con2020)',loc='right',fontsize=10)

plt.xticks([0,60,120,180,240,300,360],              [360,300,240,180,120,60,0],fontsize=14)
plt.yticks([180,150,120,90,60,30,0],              [-90,-60,-30,0,30,60,90],fontsize=14)
plt.tight_layout()

ct = ax.contour(longitude,latitude,Bg,colors='blue',levels=np.linspace(-18,0,16))
ax.clabel(ct,              inline=True,              fontsize=14,fmt='%1.1f')

ct2 = ax.contour(longitude,latitude,Bg,colors='red',levels=np.linspace(0,20,16))
ax.clabel(ct2,              inline=True,              fontsize=14,fmt='%1.1f')

divider = make_axes_locatable(ax)

#plt.savefig("JRM33-Con2020-Radial-Mag-Map.png",bbox_inches='tight')





lat_in = np.linspace(0,180,181)-90
lon_in = np.linspace(0,360,361)

latc = 0.5*(lat_in[1:] + lat_in[:-1])
lonc = 0.5*(lon_in[1:] + lon_in[:-1])


longc,latgc = np.meshgrid(lonc,latc)

Bg = Bmag(longc,latgc,model='jrm33',radial=True)



fig = plt.figure(figsize=(15,              12))
ax = plt.subplot(111)
ax.set_aspect(1.0)
ax.set_xlabel('SIII West Longitude ($^\circ$)',              fontsize=14)
ax.set_ylabel('SIII Latitude ($^\circ$)',              fontsize=14)
ax.set_xlim([360,0])        
plt.title('Iso-contours of the strength of the radial component of the field (JRM33 + Con2020)',loc='right',fontsize=10)


ct = ax.contour(longc,latgc,Bg,colors='blue',levels=np.linspace(-18,0,19))
ax.clabel(ct,              inline=True,              fontsize=14,fmt='%1.1d')

ct2 = ax.contour(longc,latgc,Bg,colors='red',levels=np.linspace(0,20,21))
ax.clabel(ct2,              inline=True,              fontsize=14,fmt='%1.1d')

divider = make_axes_locatable(ax)



lon_inZ=np.array([0,0,0,0])+85+0.5
lat_inZ=np.array([-80,-60,-40,-20])+9+0.5


longcZ,latgcZ = lon_inZ,lat_inZ

BgZ = Bmag(longcZ,latgcZ,model='jrm33',radial=True)


# plt.scatter(lon_inZ,lat_inZ)


plt.show()


# plt.imshow(Bg2d)
# plt.show()
# plt.plot(lat_in,Br2*1e-5)
# plt.plot((-1)*lat_in,Bg2d[:,360-160])






